Emergent structure in a dipolar Bose gas in a one-dimensional lattice 
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We consider an ultracold dipolar Bose gas in a one-dimensional lattice. For a sufBciently large 
lattice recoil energy, such a system becomes a series of non-overlapping Bose-Einstein condensates 
that interact via the long-range dipole-dipole interaction (ddi). We model this system via a coupled 
set of non-local Gross-Pitaevskii equations (GPEs) for lattices of both infinite and finite extent. We 
find significantly modified stability properties in the lattice due to the softening of a discrete roton- 
like mode, as well as "islands" in parameter space where biconcave densities are predicted to exist 
that only exist in the presence of the other condensates on the lattice. We solve for the elementary 
excitations of the system to check the dynamical stability of these solutions and to uncover the 
nature of their collapse. By solving a coupled set of GPEs exactly on a full numeric grid, we show 
that this emergent biconcave structure can be realized in a finite lattice with atomic ^'^Cr. 



I. INTRODUCTION 

Recent progress on the experimental realization of ul- 
tracold dipolar quantum fluids is providing an unprece- 
dented opportunity to study these systems in detail. In- 
teresting dipolar effects have been demonstrated in sys- 
tems with modest dipolar interactions, such as Bose- 
Einstein condensates (BECs) of ^^Cr [l], Addition- 
ally, the achievement of a near-quantum-degenerate gas 
of fermionic ^OR^^Rb [|, the transfer of bosonic ^iR^^Rb 
into its rovibrational ground state [3] and the cooling and 
trapping of atomic Dy ^] shows promise of a rich future 
for this field as these species host large dipole moments 
that may demonstrate new physics both in the mean-field 
regime and beyond. 

Because of the anisotropic nature of the dipole-dipole 
interaction (ddi), the physics of a dipolar system de- 
pends strongly on the geometry of the trap in which it is 
held. For example, inelastic scattering processes of both 
bosonic and fermionic species are predicted to be highly 
suppressed in tighter, quasi-two dimensional (q2D) traps 
when the trap is applied along the polarization axis of 
the dipoles |6l-[l0|. This suppression leads to more sta- 
ble, longer-lived many-body systems of reactive species. 
Additionally, it was shown using a BEC of atomic ^^Cr 
that tighter confinement in the polarization direction en- 
ergetically stabilizes a dipolar BEC (DBEC) against col- 
lapse [ll|. So, tight trapping along the polarization axis 
is necessary to obtain stable, high density dipolar quan- 
tum fluids. 

Such a trap is realizable in a one-dimensional (ID) op- 
tical lattice, where a laser is reflected onto itself and high 
and low intensity regions are formed by its interference 
pattern. The presence of the lattice brings up an interest- 
ing point regarding the physics of such a system. While 
the ddi is anisotropic, it is also long-range, scaling as 
and if the lattice spacing is sufficiently small then 
the effect of the ddi is non-negligable between the lattice 
sites. For example, interlayer superfluidity is predicted 
to exist in two adjacent layers of polar fermions and 
scattering in the 2D plane is predicted to be significantly 
modified by the presence of a weakly bound state of 



dipoles in adjacent layers [Tj]- Dramatic effects are pre- 
dicted for layers of bosons, as well, for both q2D [13, HI] 
and radially trapped [H, [13] lattice sites. The presence 
of the lattice is predicted to significantly alter the disper- 
sion via the softening of a roton-like mode in the system, 
and thus to alter the stability properties of the Bose gas. 

In this work, we consider a gas of bosonic dipoles in 
a ID lattice with the dipoles polarized along the lat- 
tice axis, so that the system is cylindrically symmetric. 
Assuming that the lattice recoil is sufficiently large, we 
model the potentials of the individual sites as cylindri- 
cally symmetric harmonic traps. At ultracold tempera- 
tures, this leads to a lattice of non-overlapping DBECs 
coupled by the long-range ddi. We study the stability 
of this system both for an infinite and finite ID lattice. 
Additionally, we find regions in parameter space where 
biconcave structure is predicted to exist that is emergent 
in the lattice system, in other words, that does not exist 
in a single condensate. To ensure the accuracy of our 
results, we calculate the elementary excitations of the 
system and use them to determine whether our solutions 
are dynamically stable. In doing so, we map the structure 
and stability of ID lattice of purely dipolar DBECs. 



II. FORMALISM 

We consider an ultracold, dilute gas of bosonic dipoles 
in a ID optical lattice in the z-direction with lattice 
spacing diat. If the lattice is sufficiently deep, it can 
be modeled by a series of iViat harmonic traps, where 
each site is described by a cylindrically symmetric po- 
tential Uj{r) = iMw^ (p2 + a2(z - jdiat)^), where M 
is the mass of the individual bosons and A = ujzjujp is 
the trap aspect ratio. This system is well-described by 
the coupled set of non-local Gross-Pitaevskii equations 
(GPEs) 
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where hj(r) is the non-interacting, or single-particle 
Hamiltonian 



operator, 



2M 



(2) 



^j(r) is the condensate wavefunction at site j, j is an 
integer and ^.j is the corresponding chemical potential. 
Without the presence of the long-range dipole-dipole in- 
teraction (ddi), these Mat equations would be indepen- 
dent. The ddi couples the equations through the mean- 
field potentials (/'^(r), given by the convolution 



0^,(r) = j dv'Va{v-v')n^{v') 



(3) 



where nj(r) — |\I'j(r)p is the density of the condensate 



occupying the j**^ site with norm J dr'nj{r') = Nj, N. 



is the condensate number for site j and Vrf(r — r'j if 
the two-body ddi potential for dipoles polarized along z 
given by [lal 



V(r - r') = d' 
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where d is the dipole moment of the bosons and ^r-r' 
is the angle between r — r' and z. A description of the 
fully-condensed, stationary state of this system of dipo- 
lar Bose-Einstein condensates (DBECs) is then given by 
the set of solutions {\E'j(r), /i^} that minimize the energy 
functional corresponding to Eq. ([T]), given by 

j 

-, Mat ^ 

Generally, a full description of a dilute EEC of inter- 
acting atoms includes contact interactions given by the 
pseudo-potential Vc(r — r') = g6{r — r') where g oc 
and flg is the s-wave scattering length of the atoms. This 
interaction is short-range, and results in the mean-field 
potential 4>i{v) = (7|^'j (r)p. Modeling a system of non- 
overlapping BECs in a ID lattice interacting only via 
contact interactions results in a set of uncoupled GPEs. 
While the interplay of contact and ddi interactions is pre- 
dicted to produce interesting effects [l^ that would likely 
be modified by the presence of the lattice, we set = 
in this work to illuminate purely dipolar effects. Because 
of its long-range nature, the ddi does not produce a sim- 
ple mean-field like the contact interaction, and requires 
particular attention. 

In practice, the dipolar mean-field is calculated in k- 
space to eliminate the problems associated with the di- 
vergence of the ddi in real-space. To do this, the Fourier 
transforms of the densities nj{r) and the ddi Vrf(r — r') 
must be calculated. Where J- is the Fourier transform 



(k) =T[nj{r)] = / drnj{r)e 



(6) 



Here, it will prove useful to define the shifted densities 
j/j(r) = nj{rj) where rj = {p,z — jdut}, so that all ^■j(r) 
are formally centered about the origin. Then, i'j{r^j) = 
nj{r), and we can write nj(k) = J^[i^j{r-j)]. With some 
simple manipulation, this expression reduces to 



(7) 



So, the fc-space density of the DBEC at site j can be 
rewritten as the Fourier transform of njij) translated 
into the local set of coordinates, with an additional expo- 
nential term accounting for this spatial translation. Now, 
by the convolution theorem, the mean-field contribution 
from the DBEC at site j can be written as 

0^,(r)=^-i [y,(k)n,(k)] (8) 
where Vd(k) is the Fourier transform of the ddi [20| . 
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and add — is the characteristic dipole length. In 

this work, we calculate Vj{v) directly by calculating 
the shifted condensate wavefunctions 3>j(r) such that 
i'j(r) — |$j(r)p and account for the spatial separation 
of the DBECs, or the presence of the lattice, with the 
expression given in Eq. ([T]). So, the wavefunctions vE'j(r) 
and $j(r) are related by ^j{v) = '^j{Yj). 



III. WAVEFUNCTION ANSATZ 

For a single DBEC, calculating the mean-field en- 
ergy on a full numeric grid has proven fruitful [2ll [2^ , 
however, this method is very computationally expensive 
when considering multiple interacting DBECs, both in 
real- and fc-space. In real-space, the convolution integral 
for the dipole-dipole mean-field must be done directly, 
where there is no divergence if the condensates do 
not overlap. In fc-space, the grid must be large enough to 
resolve the entire lattice because of the e*'^^ depen- 
dence of the fc-space densities. To avoid these problems, 
we consider solutions of the form $j(r) = ^j{p)xj{'^) 
where 
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FIG. 1: (color online) The differences in energy of a DBEC 
in a trap with aspect ratio A = 7 as a function of interaction 
strength between that calculated using the ansatz given in 
Eq. (|10|l and that calculated exactly on a full numeric grid. 
The blue dotted line shows the energy difference calculated 
using no 2""* order h.o. wavefunction and a fixed axial width 
Iz ~ o-z, the red dashed line shows the energy difference using 
the same wavefunction but with Iz treated variational and the 
black solid line shows the energy difference including the 2'^'* 
order h.o. wavefunction where the relative amplitude A2 and 
Iz are treated variationally. 
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where H2 (x) — Ax^ — 2 is the second Hermite polynomial. 
This ansatz includes the zeroth and second harmonic os- 
cillator wavefunctions with variable width and relative 
amplitude. Plugging this ansatz into the GPE and inte- 
grating out the z-dependence results in a modified GPE 
in the radial coordinate p that also depends on the widths 
Iz^j and the relative amplitudes A2^j of the axial wave- 
functions, but not the z-coordinate explicitly. We derive 
this modified GPE for a single DBEC, given by Eq. (|X2|) . 
in Appendix \K\ 

To test the ansatz given in Eq. (ITUl) . we apply it to 
the well known system of a single DBEC in a harmonic 
trap. This system was predicted to exhibit, for certain 
trap geometries and ddi strengths, biconcave structures 
where the maximum density of the DBEC exists not in 
the center of the trap, but in a ring about the center of 
the trap Q . For example, such structure is predicted to 
exist in a trap with aspect ratio A = 7 for ddi strengths 
near the stability threshold. 

Figure [T] compares the total energies of a single DBEC 
in a trap with A = 7 as a function of the ddi strength 
{N — l)gd/ap, where Op — y'h/Mujp is the radial har- 
monic oscillator length and cjd = ^^^^ is the ddi 
coupling, for various restrictions placed on the variational 
parameters of the axial wavefunction. Plotted is the en- 
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FIG. 2: (color online) The values of the axial wavefunction 
parameters that, together with the radial wavefunction cal- 
culated on a grid, minimize the energy of a single DBEC in 
a trap with aspect ratio A = 7. The blue dotted line shows 
the result for Iz — az, the red dotted line shows the values of 
Iz when it is treated variationally and A2 = and the black 
solid line shows the values of Iz when it and A2 , the black dot- 
ted line (marked by the right vertical axis) are both treated 
variationally. 



ergy difference (i? — i?fuu)/i?fuii, where Ei^n is the energy 
calculated by solving the GPE exactly (within strict nu- 
merical precision) on a full numeric grid in p and z. The 
blue dotted line shows the energy of the DBEC when 
A2 = and Iz is fixed to be the axial harmonic oscillator 
length, ttz = \/ h /Mujz , the red dashed line shows the 
energy when A2 = and Iz is treated variationally, and 
the black line shows the energy when A2 and Iz are both 
treated variationally. Clearly, the full variational treat- 
ment is much more accurate than the cases where the 
second harmonic oscillator wavefunction is not included 
{A2 = 0). Indeed, it stays within 1% of the exact en- 
ergy for all values of {N — l)gd/ap for which the DBEC 
is stable. We find this to hold true for larger trap aspect 
ratios, as well. Figure [5] shows the values of the varia- 
tional parameters for the same cases as in figure [T] In 
this figure, the left vertical axis labels Iz/az and the left 
vertical axis labels A2, shown by the black dots. 

Beyond energetics, this ansatz also predicts semi- 
quantitatively the structure and stability of a single 
DBEC. An interesting feature of the biconcave struc- 
ture predicted in Ref. [1^ is that it exists in "islands" 
of parameter space, defined by {N — l)gd/ap and A. Fig- 
ure [3] shows this structure/stability diagram for a single 
DBEC calculated using a) A2 — and Iz/ciz = 1, b) A2 
and Iz variational and c) a full numeric grid in p and z. 
Interestingly, the biconcave islands are present in each 
diagram and occur for almost exactly the same values of 
{N — \)gd/ap. They are, however, shifted in A, moving 
to smaller values as more restrictions are placed on the 
condensate wavefunction. The diagram for the full vari- 
ational ansatz (b) qualitatively matches that of the full 
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FIG. 3: (color online) Structure/stability diagram for a sin- 
gle DBEC. The colored regions indicate a dynamically stable 
condensate, and the pink (darker) regions indicate parameters 
for which the DBEC has biconcave density, (a) and (b) are 
calculated using the ansatz for the axial wavefunction given 
in Eq. (|10p and (c) is calculated using a full numeric grid. For 
(a), A2 = and Iz/a^ = 1, and for (b), A2 and h are treated 
variationally. 

numeric grid. Thus, we expect that this ansatz will give 
physically meaningful results, if not quite quantitative 
ones. 

As we discuss below, a key benefit of this ansatz for 
the ID lattice system is that it is analytic in z. An- 
other such ansatz that has this property is that of cor- 
related Gaussians, which have been shown to reproduce 
the results of full numeric calculations for DBECs quite 
well . However, we applied this ansatz to the lattice 
system and found that it is numerically unstable with the 
minimization techniques used here. 



to consider is that of an infinite lattice, with A'lat — > 00. 
This approximation introduces a discrete invariance to 
the system so that we can set <&j(r) — ^j'(r) for all 
Thus, we can neglect, for the time being, the indexing of 
the wavefunctions and let $j(r) — > <i>(r) for all j. Then, 
the mean field potential at any site is, from Eq. ([7]), given 
by 

CO „ 

Mr)^ E / «(k)n(k)e"=^'^'-^e-* '- (13) 

j = -oo 

where 

e i'^i'i 

n(k) = np{kp) ^^jq 

X (1 + (^A, - + ^2) kill + lA.ktlt^'^ .(14) 

We can manipulate the infinite sum in Eq. (|13p to 
give [H 

J — — 00 j — ~oo 

The term that accounts for the infinite lattice can there- 
fore be written as a Dirac comb in k^ with spacing 27r/(iiat 
between peaks. Inserting this expression into Eq. (fT3)) 
gives the mean-field potential 

Mr) = 2-^J'2B [Finf{kp)Mkp)] , (16) 

«lat 

where 2gdFi-a_i{kp) / dut is the effective fc-space ddi for the 
infinite lattice and i^inf (^p) is given by 



IV. INFINITE LATTICE 



With confidence in the ansatz given in Eq. (jlOp . we now 
apply it to the ID lattice system. An interesting example 
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The GPE for an infinite lattice of interacting DBECs 
is reduced to a single GPE in the radial coordinate p 
where all of the axial dependence of the wavefunction is 
captured by the variational parameters A2 and Iz. 

We study the structure and stability of this infinite 
lattice of interacting DBECs by solving the modified 



GPE for the system (applying conjugate gradients [2^ 
to minimize the corresponding energy functional) and 
studying the Bogoliubov de Gennes (BdG) excitations. 
We find that the sum in Eq. ([T7| is sufficiently con- 
verged if a cutoff jcut is applied to the index j such that 

jcut » dlat/27r/2. 
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FIG. 4: (color online) Structure/stability diagram for an in- 
finite lattice of DBECs in traps with aspect ratio A = 10 as 
a function of lattice spacing di^t / a,z and interaction strength 
{N — l)gd/ap. The colored region indicates dynamic stability, 
while the pink (darker) regions indicate parameters where the 
DBECs have biconcave density. The inset shows an isodensity 
plot of a DBEC with biconcave density. 



FIG. 5: (color online) Structure/stability diagram for an in- 
finite lattice of DBECs in traps with aspect ratio A = 20 as 
a function of lattice spacing diat/oz and interaction strength 
(A^ — l)gd/ap. The inset shows a close-up of the diagram at 
the parameters indicated. The pink (darker) region in the 
inset indicated parameters where the DBECs have biconcave 
density. An isodensity plot of a DBEC with biconcave density 
is shown in this inset. 



Consistent with other results [T3, [13 , we find that the 
presence of the lattice serves to destabilize the system 
due to the softening of a discrete roton-like mode in the 
system. For a single DBEC in a trap, tight axial confine- 
ment aligns the dipoles so that they are predominately 
repulsive and, for sufficiently low densities or interactions 
strengths, stabilizes the condensate. In the presence of a 
ID lattice, the attraction from the dipoles at other lat- 
tice sites extends the condensate in the axial direction, 
increasing the integrated axial density and, ultimately, 
making the system less stable. This destabilization is 
made less dramatic as diat is increased. 

To study the structure and stability of the infinite lat- 
tice, we choose specific trap aspect ratios and explore the 
parameter space defined by {N — l)gd/ap and diat. Fig- 
ure |4] shows the region of dynamic stability for an infinite 
lattice of DBECs in traps with A — 10. For lattice spac- 
ings d\s±/az ^ 5, the condensate wavefunctions at adja- 
cent sites overlap and the strong dipole-dipole attraction 
leads to complete instability. In this figure, the colored 
regions indicate dynamic stability and the pink (dark) 
regions indicate parameters at which the DBECs exhibit 
biconcave density. As diat/oz is increased, the diagram 
approaches that given by a line at A = 10 in figure [3] for a 
single DBEC. However, for smaller lattice spacings, a sec- 
ond biconcave island appears. Without the presence of 
the lattice, biconcave structure would not exist for these 
parameters. Thus, this structure is "emergent" in the 
lattice system. The inset in figure |¥] shows an isodensity 
plot of a DBEC with biconcave density. 

Figure [5] shows the region of dynamic stability up to 
lattice spacings of diat/oz = 80 for a infinite lattice of 
DBECs in harmonic traps with A = 20. Here, the con- 
vergence of the stability line to [N — l)gd/ap ^ 550 is 
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FIG. 6: (color online) Stability lines for an infinite lattice of 
DBECs for aspect ratios A = 50, 100, 150 as a function of lat- 
tice spacing diat/iz and interaction strength (A'^ — l)gd/ap- 
The parameters beneath the lines are dynamically stable, 
while those above the lines are dynamically unstable. 



clear. The inset shows a close-up view of the diagram 
where a biconcave island is predicted to exist. As the as- 
pect ratio is increased, the values of interaction strength 
{N — Vjgd/cLp that the biconcave islands span becomes 
relatively smaller compared to the asymptotic value of 
the stability line. Figure [5] shows the stability lines for 
infinite lattices with aspect ratios A = 50, 100, 150. We 
find stability islands that exist at the stability thresh- 
old within the lattice spacings rfiat/a^ = 6 to 10 for all 
of these aspect ratios. Because they are so narrow in 
{N — l)gd/ap, though, they are not included in this plot. 



6 



By working in the p- and z-coordinates, a cylindrical 
symmetry is assumed. However, it was shown in (23| 
that DBECs with biconcave densities are dynamically 
unstable to angular modes, or quasiparticles. While the 
method used here is sensitive to dynamic instabilities 
that are purely radial, an extra step must be taken to 
detect angular instabilities. 



Bogoliubov de Gennes Equations 

The Bogoliubov de Gennes (BdG) equations describe 
the low-lying quasiparticles of the condensate. They are 
derived by inserting the ansatz 



-ifit 



'(18) 

into the time-dependent GPE and linearizing about 5, 
assuming that 5-^1. Here, m is the quantum num- 
ber describing the projection of angular momentum of 
the quasiparticle onto the z-axis. In general, the energy 
eigenvalues of the BdG modes {u, v*} can be written as 
Lo — LOji + iioi, where ujn and w/ are purely real. When 
all Loi = 0, the system is dynamically stable. However, 
when some wj ^ 0, the system is dynamically unstable, 
and the quasiparticle amplitude grows exponentially in 
time on a time scale ~ l/ujj. 

Like the single DBEC in we find that the bicon- 
cave structures in the infinite lattice are, for some critical 
density or ddi, dynamically unstable to angular quasi- 
particles with TO > 2. Ref. [21] shows that this angular 
instability leads to angular collapse, or collapse with an- 
gular nodes, of the biconcave DBECs. A measurement 
of the character of collapse, whether it be radial or angu- 
lar, then provides a tool to map the structure along the 
stability threshold of the system. 

In our analysis, we found exotic ground state densities 
very close to the stability threshold, like those found for 
a finite lattice in Ref. jl^]. These solutions host multi- 
ple radial density oscillations, however, we find they are 
dynamically unstable and are thus unlikely to be exper- 
imentally observable. 



V. FINITE LATTICE 

While the infinite lattice of DBECs provides a clear, 
simple example of emergent structure in this system, it is 
a difficult system to realize experimentally. In a realistic 
experiment, the lattice has a finite extent and the occu- 
pations of the sites vary from site to site. To model this 
more realistic lattice system, we consider an odd num- 
ber of occupied lattice sites indexed by j G [— jiat,Jiat] 
where j'lat — (Mat — l)/2 with particle number given 
by a Gaussian distribution, Nj — A^maxGxp [— (j'/j'iat)^], 
where A^max is the particle number in the condensate in 
the center of the lattice at site j = 0, a.nd the outer-most 
sites have particle number A^max/e |16| . 
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FIG. 7: (color online) Radial densities at 2: = of a DBEC 
with (A'^ — l)gd/ap — 550 in a trap with aspect ratio A = 50. 
The blue dashed line shows the density of a DBEC in a single 
harmonic trap, and the red solid line shows the density of a 
DBEC in the center site of a ID lattice with nine occupied 
sites (jiat = 4). This DBEC exhibits biconcave structure, 
while the single DBEC does not, demonstrating the emer- 
gence of this structure in the lattice system. These densities 
were calculated by solving the GPE (coupled GPEs) exactly 
on a full numeric grid. 



Instead of using an analytic form for the axial parts of 
the condensate wavefunctions, we solve the coupled set 
of GPEs given by Eq. ([T|) on a full grid (large enough to 
encapsulate the entire lattice) in p and z for each con- 
densate. We find good convergence by using the conju- 
gate gradients method to minimize the full energy func- 
tional of the system 26] . Additionally, to ensure numer- 
ical precision we apply a cutoff to the ddi in p and z 
so that a relatively small grid can be used while elimi- 
nating the effects of artificial "image" condensates that 
are present due to the use of the FFT algorithm in our 
calculation | 19j ]. 

As an example, we consider a lattice with trap as- 
pect ratios A = 50, jiat = 4 (corresponding to 9 oc- 
cupied lattice sites), lattice spacing diat = 8a^ and 
(Amax — ^)9d/ap — 550 on a numeric grid of size 
[A^z,A^p] = [1024,128]. Figure [3 shows the density at 
z = of a DBEC at the center of the lattice (j — 0) and, 
for comparison, the density of a DBEC with the same 
trap aspect ratio and ddi strength (A^ — l)gd/ap — 550 
but without the presence of the lattice. While the DBEC 
in the single trap does not exhibit biconcave structure, 
the DBEC in the lattice does, showing that this emergent 
structure in the lattice system is present not only in the 
infinite lattice system, but also in the experimentally re- 
alistic system of a finite lattice with variable occupancy. 
Indeed, such a system is realizable with atomic ^^Cr, hav- 
ing a permanent magnetic dipole moment of ^ = Qps 
where pB is the Bohr magneton, axial harmonic oscil- 
lator frequencies of = 27r x 30 kHz and a maximum 
condensate occupancy of A'max — 77 x 10'^ atoms. 
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VI. CONCLUSION 

In conclusion, we have mapped the structure and sta- 
bihty of a lattice of interacting, purely dipolar DBECs. 
By asserting an analytic form for the axial part of the 
condensate wavefunctions (Eqs. (ITm - (fT2]) '). we derive a 
simple, modified GPE for the radial part of the wave- 
functions when the lattice is infinite. We find isolated 
regions ( "islands" ) in the parameter space defined by the 
lattice spacing and the ddi strength where the DBECs are 
predicted to exhibit biconcave densities, where the max- 
imum density exists not in the center of the trap but in 
a ring about the center of the trap. To model a more ex- 
perimentally realistic system, we consider a finite lattice 
with varying condensate number and solve the coupled 
set of GPEs exactly on a full numeric grid. In doing so, 
we show that this emergent biconcave structure should 
be observable in a finite ID lattice of DBECs of atomic 
52Cr. 



Acknowledgments 

We acknowledge the financial support of the NSF and 
the DOE, and would like to thank Juliette Billy and 
Jonas Metz for their enlightening discussions. 



Appendix A: Modified GPE using 0* and 2""^ 
harmonic oscillator wavefunctions 



Consider the ansatz given by Eqs. (fT0 |) - (fT2|) for a sys- 
tem with a single harmonically trapped DBEC, so the 
indexing of the condensate wavefunction can be ignored 
and we can simply write $(r) — ip(p)x{z). We derive the 
modified GPE by multiplying the (dimensionless) GPE, 



1 



V" + C/(r) + Mr) - } *(r) = 0, 



(Al) 



by xi^) and integrating over z. This operation gives the 
modified GPE, 



9d_ 

^2 



71p(fcp)Fcff ( ^ 



■ V'(p) = 0, 

(A2) 

where /loff(p) is the effective single-particle Hamiltonian, 



l-v2 I 1, ,2 „2 



/icfr(p) = --V^ + -cjy-M 



1 



1 



1 



2,2 



and i^cfi(a;) is given by 



,(A3) 



Fesix) = 7 I 1 + V2A2 (3x^ - 1) + -Al (3 + 5x^ + 6x'^) + -^A\ (l + <dx^ + 42a;'' \2x^\ 

{\ + Alf \ V ; 4 2V ^ 4V2 ^ ' 

^ ^A\ (41 + 3x2 ^ ^34^2 ^ go^4 ^ 8^6)) _h/lu^A2 (2 -f A^^ 
54 2 

2 (^V2 + 2A2) (1 + A^) + (3 + 4\/2A2 + 5^2) + |^x/2 + 2A2) a;^ + i A^x^ 



xe^ Erfc [x] (A4) 



and Erfc[x] is the complimentary error function. The 
corresponding mean-field energy due to the ddi is then 



2D 



(A5) 
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